function [err] = EPSC_mini(n)
%
% EPSC_mini: mini analysis
% based on trace_view: view individual or selected traces, with zooms etc.
% 8/13/01 Paul B. Manis.
% 6/20/02 Updated to include the Clements and Bekkers method
% and offer a choice using a gui.
%
% 6/7/2007. Added ability to select a time window within the trace for
% analysis.
% Added a max amplitude criterion for the template analysis (if data in
% selected window exceeds the max amplitude, + or -, the event is
% rejected). This allows rejection of noise spikes and stimulated EPSCs.
%

global STACK
decay_start = 0.85; % position for decay fitting...
fit_re_peak = 1; 	% do exp fits from peak of trace rather than "0" time
persistent im

err = 0;



switch(nargin)
    case 0   % handle the display...
        % do parameter entry ...
        
        if(exist('mini_par.mat', 'file') == 0) % read the file... if it exists
            % Otherwise initialize the structure and save it
            % for standard thresholding analysis.
            par.posneg = 1;	% 1 is negative, 2 is positive.
            par.posstd = 3;
            par.negstd = 3;
            par.maxamp = 500;
            par.tpre = 2;
            par.tpost = 5;
            par.nbins = 32; % histogram number of bins.
            par.lpf = 2000; % lpf filter
            par.channel = 2;
            % and for template analysis
            par.t_rise = 1;
            par.t_decay = 5;
            par.t_thresh = 3.0;
            par.t_sign = 1;
            par.t_lpf = 2000;
            par.t_nbins=32;
            par.tstart = 0; % analsis window start/end times
            par.tend = 10000;
            par.t_maxamp = 300; % template max amplitude - rejection criteria
            par.template = 1;
            par.ntau = 4;
            save('mini_par.mat', 'par');
        else
            load('mini_par.mat', 'par');
            if(~isfield(par, 'channel')) %#ok<NODEF>
                par.channel = 2;
            end;
            if(~isfield(par, 'disppre'))
                par.disppre = 5;
                par.disppost = 40;
            end;
            if(~isfield(par, 'tstart'))
                par.tstart = 0;
                par.tend = 10000;
            end;
            if(~isfield(par, 't_maxamp'))
                par.t_maxamp = 500;
            end;
            if(~isfield(par, 'template'))
                par.template = 1;
            end;
            if(~isfield(par, 'ntau'))
                par.ntau = 4;
            end;
        end;

        % ok, minipar is set up.

        % now set up the gui window
        hm = findobj('tag', 'EPSCAnalysis');
        if(isempty(hm))
            open('minifig.fig');
            hm = findobj('tag', 'EPSCAnalysis');
            datac('addwindow', 'EPSCAnalysis');
        end;
        figure(hm);

        set_window_arg('Mini_channel', par.channel);
        set_window_arg('Mini_disppre', par.disppre);
        set_window_arg('Mini_disppost', par.disppost);

        set_window_arg('Mini_std_pre', par.tpre);
        set_window_arg('Mini_std_post', par.tpost);
        set_window_arg('Mini_std_amp', par.maxamp);
        set_window_arg('Mini_std_stdp', par.posstd);
        set_window_arg('Mini_std_stdn', par.negstd);
        set_window_arg('Mini_std_bins', par.nbins);
        set_window_arg('Mini_std_lpf', par.lpf);
        set_window_arg('Mini_std_sign', par.posneg, 1);

        set_window_arg('Mini_temp_rise', par.t_rise);
        set_window_arg('Mini_temp_decay', par.t_decay);
        set_window_arg('Mini_temp_maxamp', par.t_maxamp);
        set_window_arg('Mini_temp_thresh', par.t_thresh);
        set_window_arg('Mini_temp_bins', par.t_nbins);
        set_window_arg('Mini_std_lpf', par.t_lpf);
        set_window_arg('Mini_temp_sign', par.t_sign, 1);

        set_window_arg('Mini_tstart', par.tstart);
        set_window_arg('Mini_tend', par.tend);
        set_window_arg('Mini_template', par.template, 1);
        set_window_arg('Mini_ntau', par.ntau);
        return;

    case 1 % we have an input argument... must be call from the miniEPSC windowlysis..
        % so read the window
        load('mini_par.mat', 'par');
        if(~isfield(par, 'channel')) %#ok<NODEF>
            par.channel = 2;
        end;
        if(~isfield(par, 'disppre'))
            par.disppre = 5;
            par.disppost = 40;
        end;
        if~(isfield(par, 'tstart'))
            par.tstart = 70;
            par.tend = 90;
        end;
        if(~isfield(par, 't_maxamp'))
            par.t_maxamp = 500;
        end;
        if(~isfield(par, 'template'))
            par.template = 1;
        end;
        if(~isfield(par, 'ntau'))
            par.ntau = 1;
        end;
        par.channel = get_window_arg('Mini_channel', par.channel); % always get these...
        par.disppre = get_window_arg('Mini_disppre', par.disppre);
        par.disppost = get_window_arg('Mini_disppost', par.disppost);
        par.tstart = get_window_arg('Mini_tstart', par.tstart);
        par.tend = get_window_arg('Mini_tend', par.tend);
        par.template = get_window_arg('mini_template', par.template, 1);
        par.ntau = get_window_arg('Mini_ntau', par.ntau);
        if par.template == 5
            STACK{1} = im;
        end;
        switch(n) % find out what the command is.
            case 1 % do the standard analysis - first read the parameters and write to the par file

                par.tpre = get_window_arg('Mini_std_pre', par.tpre);
                par.tpost = get_window_arg('Mini_std_post', par.tpost);
                par.maxamp = get_window_arg('Mini_std_amp', par.maxamp);
                par.posstd = get_window_arg('Mini_std_stdp', par.posstd);
                par.negstd = get_window_arg('Mini_std_stdn', par.negstd);
                par.nbins = get_window_arg('Mini_std_bins', par.nbins);
                par.lpf = get_window_arg('Mini_std_lpf', par.lpf);
                par.posneg = get_window_arg('Mini_std_sign', par.posneg, 1); % 0 is negative, 1 is positive.
                par.t_disp = get_window_arg('Mini_std_display', 0, 1);

                save('mini_par.mat', 'par'); % update the file
                % now do the analysis

                sf = get_window_arg('CTL', 1, 1);
                for i = 1:length(sf)
                    [err, im] = find_transients(1, sf(i), decay_start, fit_re_peak); % direct call...
                end;
                return;

            case 2
                par.t_rise = get_window_arg('Mini_temp_rise', par.t_rise);
                par.t_decay = get_window_arg('Mini_temp_decay', par.t_decay);
                par.t_thresh = get_window_arg('Mini_temp_thresh', par.t_thresh);
                par.t_nbins = get_window_arg('Mini_temp_bins', par.t_nbins);
                par.t_lpf = get_window_arg('Mini_temp_lpf', par.t_lpf);
                par.t_sign = get_window_arg('Mini_temp_sign', par.t_sign, 1);
                par.t_disp = get_window_arg('Mini_temp_display', 0, 1);
                par.t_maxamp = get_window_arg('Mini_temp_maxamp', par.t_maxamp);
                save('mini_par.mat', 'par'); % update the file
                % now do the analysis

                sf = get_window_arg('CTL', 1, 1);
                for i = 1:length(sf)
                    [err, im] = find_transients(2, sf(i), decay_start, fit_re_peak); % direct call...
                end;
            otherwise
        end;
        
    otherwise
end;

%----------------------------------------------------------------------------------------

function [err, im] = find_transients(cmd, sf, decay_start, fit_re_peak)
% main analysis routine. Finds transients and does some basic
% analysis assuming that they are supposed to be EPSCs.
% 6/20/02 - modified to use a real gui window
% % cmd determines analysis method
% if 1, then use standard thresholding
% if 2, then use Clements and Bekkers method of template matching
%
global DFILE ALLCH
global FILTERED_TRACE TRACE
global CONTROL
persistent analyzeflag
err = 0;
if(~exist('analyzeflag', 'var'))
    analyzeflag = 0;
end;

if(analyzeflag == 1)
    analyzeflag = 0;
    return;
else
    analyzeflag = 1;
end;


grpos = zeros(9,4);
grcxpos = grpos;
for x = 1:3
    for y = 1:3
        lx = 0.05+(x-1)/3;
        by = 0.05+(2/3)-(y-1)/3;
        xw = 0.27;
        yw = 0.27;
        rhs = lx + xw;
        top = by + yw;
        cx=0.1;
        cy=0.07;
        grpos(x+(y-1)*3,:) = [lx, by, xw, yw];
        grcxpos(x+(y-1)*3,:) =[rhs-cx, top-cy, cx, cy];
    end;
end;
%fprintf(1, 'findtransients 1');

if(nargin < 1 || nargin > 4)
    QueMessage('EPSC_mini/find_transients: check arguments in call\n');
    analyzeflag = 0;
    return;
end;

load('mini_par.mat', 'par'); % get the parameters (passed via disk.. unusual but safe)

% set up the common values between the protocols.

switch(cmd)
    case 1
        SIGN = par.posneg-1;
        NBINS = par.nbins;
        LPFFILTER = par.lpf;
        iminus = par.tpre;
        iplus = par.tpost;
    case 2
        SIGN = par.t_sign-1;
        NBINS = par.t_nbins;
        LPFFILTER = par.t_lpf;
        if(par.template < 4)
            iminus = par.t_rise*3;
            if(iminus < 0.5)
                iminus = 0.5;
            end;
            ntau = par.ntau;
            if(par.t_decay*ntau < 2) % less than 3 msec is a problem though
                ntau = 2/par.t_decay; % make it a minimum of 2 msec...
            end;
            iplus = par.t_decay * ntau;
        else
            %       ntau = 1;
            iminus = 0.5;
            iplus = par.t_decay;
        end;
    otherwise
        QueMessage('EPSC_mini/findevents: bad command');
        analyzeflag = 0;
        return;
end;


% verify file already in memory or else read it in
[DFILE] = synch_file(DFILE, sf);
if(isempty(DFILE) || (isempty(ALLCH)))
    QueMessage('EPSC_mini/find_transients: could not get DFILE', 1);
    analyzeflag = 0;
    return;
end;

if(~check_field(CONTROL(sf), 'I_Unit'))
    CONTROL(sf).I_Unit='pA';
end

if(~check_field(CONTROL(sf), 'V_Unit'))
    CONTROL(sf).V_Unit='mV';
end

if(~check_field(CONTROL(sf), 'T_Unit'))
    CONTROL(sf).T_Unit='ms';
end
if(~check_field(CONTROL(sf), 'sp_deadwin'))
    CONTROL(sf).sp_deadwin = CONTROL(sf).deadwin;
end

% get the protocol
protocol=deblank(lower(CONTROL(sf).protocol));

% determine rates for each trace.
rate = DFILE.rate.*DFILE.nr_channel/1000;
records = length(DFILE.record);
pts=DFILE.nr_points;

%if(~isempty(CONTROL(sf).psp_reclist))
%	[RL, err] = record_parse(CONTROL(sf).psp_reclist);
%end;
[RL, err] = record_parse(CONTROL(sf).reclist);
frec=min(RL);
lrec=max(RL);

% compute the time base for plotting (time is [rec, npts] array for EACH record)
k = find(DFILE.rate < 1);
DFILE.rate(k) = 50;
RATES = (DFILE.rate .* DFILE.nr_channel)/ 1000; % array of sampling rates, convert to msec


% always compute the time of each stimulus, in seconds from the start of the data
if(DFILE.mode >= 5)
    wz=DFILE.ztime;
    w=find(diff(wz) < -12*60*60); % correct for possibility that someone actually records the data across midnight rollover.... (yes, it happened. 5/16/01 with Huijie's data.)
    if(~isempty(w))
        wz(w+1:end)=wz(w+1:end)+24*60*60;
    end;
    zt = (wz-wz(1))/(60);
    cond_baseline = 5; % 5 min baseline
else
    zt = (DFILE.ztime-DFILE.ztime(1))/(60*1000);
    cond_baseline = 3;
end;
TM=zt;


% now build the window
hmini = newfigure('EPSC_minitraces', DFILE.filename);


fsize = 7;
msize = 3;
msg_pos1 = [0.85 0.15 0.15 0.07];
msg_pos2 = [0.85 0.05 0.15 0.07];


npts=DFILE.nr_points;
% time=zeros(1,npts);
k=1;
if(length(DFILE.nr_channel) > 1)
    j = k;
else
    j = 1;
end;
p=(DFILE.rate(k)*DFILE.nr_channel(j)/1000);
time=0:p:(npts-1)*p';
tmax = max(time);

QueMessage('EPSC mini analysis starting',1);

% filter the DTRACE if necessary
% Smooth the DTRACE out a bit
TRACE = ALLCH{par.channel};
TRACE = FP_artsupp(TRACE, DFILE, sf, 1); % suppress the artifacts if stimulating..... (1 is an "auto flag")


if(isempty(FILTERED_TRACE) || FILTERED_TRACE == 0)
    for i=1:records
        fsamp = 1000/RATES(i); % get sampling frequency
        fco = LPFFILTER;		% cutoff frequency in Hz
        wco = fco/(fsamp/2); % wco of 1 is for half of the sample rate, so set it like this...
        if(wco < 1) % if wco is > 1 then this is not a filter!
            [b, a] = butter(8, wco);
            TRACE(i,:) = filter(b, a, TRACE(i,:)); % filter each trace...
            if(isempty(FILTERED_TRACE) || FILTERED_TRACE == 0)
                QueMessage('Filtering');
                FILTERED_TRACE = 1;
            end;
        end;
    end;
end;

% clip the trace to a particular time window.
if(par.tend > 0 && par.tend > par.tstart)
    [a, jstart] = find(time >= par.tstart);
    jstart = min(jstart);
    [a, jend] = find(time <= par.tend);
    jend = max(jend);
    told = TRACE;
    TRACE=zeros(records, jend-jstart+1);
    for i = 1:records
        TRACE(i,:) = told(i,jstart:jend); % just cut it out...
    end;
end;

% the following is the event detection loop, where we examine each mini/event

N = size(TRACE);
ihist = NaN*zeros(1, 10000); % peak TRACE histogram data
qhist = NaN*zeros(1,10000); % total charge histogram data
isihist = ihist; % inter-event ('spike') histogram data
itau = ihist;
itaum = ihist;
iamp = ihist;
icur = zeros(1,10); % TRACE in acceptable traces
m = 0; % count acceptable traces
twin = [0 0];
im = [0 0];
stim_time = 5;
stim_list=CONTROL(sf).stim_time; % get the array
if(ischar(stim_list))
    stim_time= number_arg(stim_list);
end;
tbs = floor((stim_time(1)./RATES)+0.5)+1;
ieventtime = cell(1, records);
ieventamp = cell(1, records);
ucmd = cmd;


% get the template, because we will use it to determine the exp decay fit window (!)
[tmpl, predel] = cb_template(par.template, p, par.t_rise, par.t_decay, 0, par.ntau); % LPFFILTER); % compute base template
[temp_max, itemp_max] = max(tmpl);
i95 = decay_start*temp_max;
jmx = floor((length(tmpl) - itemp_max)*0.9) + itemp_max;
[i95f p95f] = findnear(tmpl(itemp_max:jmx), i95); % find point 95% way down falling slope
if(length(p95f) >= 2)
    p95f = p95f(2)+itemp_max - 1;
else
    p95f = p95f + itemp_max - 1;
end;
fitmeth = 'c';

for k = 1:records
    % use choice of methods
    QueMessage(sprintf('EPSC mini analysis: record %d of %d', k, records));
    %    TRACE(k, (1:(tbs-1))) = detrend(TRACE(k, (1:(tbs-1))));
    imean = mean(TRACE(k,(1:(tbs-1))));
    istd = std(TRACE(k,(1:(tbs-1)))); % mean_std('TRACE', k);
    p=(DFILE.rate(k)*DFILE.nr_channel(k)/1000);
    istdp = (par.posstd*istd); % positive and negative detection/acceptance limits.
    istdn = (par.negstd*istd);

    switch(ucmd)
        case {1} % just use straight threshold detection
            switch(SIGN)
                case 1
                    [d1, k2] = find((TRACE(k,:)-imean) > istdp); % find points > threshold
                    icand = k2(find(diff(k2) > 1)+1); % and store the resulting array
                case 0
                    [d1, k2] = find((TRACE(k,:)-imean) < -istdn); % find points > threshold
                    icand = k2(find(diff(k2) > 1)+1); % and store the resulting array
                otherwise
                    icand = [];
            end;
        case {3} % use modified with derivative...
            switch(SIGN)
                case 1
                    [d1, k2] = find((diff(TRACE(k,:))) > istdp); % find points > threshold
                    icand = k2(find(diff(k2) > 1)+1); % and store the resulting array
                case 0
                    [d1, k2] = find(diff((TRACE(k,:))) < -istdn); % find points > threshold
                    icand = k2(find(diff(k2) > 1)+1); % and store the resulting array
                otherwise
                    icand = [];
            end;
        case 2 % use template method for detection
            [icand, iscamp, ioff, crit] = ClementsBekkers(TRACE(k,:)-imean, ...
                p, par.t_rise, par.t_decay, par.t_thresh, SIGN, par.t_disp, 0, par.template, par.ntau, hmini); % LPFFILTER, 2);

        otherwise
    end;
    jl = 0; % index to last real event
    evcount = 0; % total count of events in this trace

    lj = 0;
    ieventtime{k} = NaN*ones(1, 10000);
    ieventamp{k} = NaN*ones(1, 10000);
    QueMessage(sprintf('EPSC mini analysis: %d candidate events in %d traces', length(icand), length(records)));
    %   fprintf(1, 'cand events: %d\n', length(icand));
    ncand = length(icand);
    for i = 1:ncand % go through all the candidate points
        if(analyzeflag == 0)
            return;
        end;
        j = icand(i);
        % check for overlap on the falling phase of the events
        poverlap = 1;
        if(i < ncand) % check cases for overlap on falling side
            if(j + floor(iplus/p) < icand(i+1)) % no overlap with next putative event
                poverlap = 0;
            end;
        else % last candidate event - see that it fits in the trace duration
            if(size(TRACE,2) > icand(i)+floor(iplus/p)) % and fits inside the trace
                poverlap = 0;
            end;
        end;

        noverlap = ((j - jl) < floor(iminus/p));% also make sure about the rising phase.

        if (~noverlap  && ~poverlap )
            prewin = j - (5:(floor(iminus/p)));
            twin = (-iminus:p:iplus); % generates a time base for this plot
            iwin = (1:length(twin)) + (j-floor(iminus/p));
            % criteria for acceptance:
            %    must fit in the data window completely,
            %    there must be no large events greater than maxamp
            if(min(prewin) > 0 && ...
                    min(iwin) > 0 && ...
                    max(iwin) < npts) % & ...
                %               all(abs(TRACE(k,iwin)-imean) < par.maxamp)); % & ...
                %             fprintf(1, '%d fits in window\n', i);
                twin = (-iminus:p:iplus); % generates a time base for this plot
                qwin = j + (-1:floor(2/p));
                iwin = (1:length(twin)) + (j-floor(iminus/p));
                bl=mean(TRACE(k,prewin));
                iplot = TRACE(k,iwin)-bl; % remove the baseline once and for all
                %iplot = detrend(TRACE(k, iwin)); % try tihs instead...
                go = 1; % final inclusion criterion - no junk going the other way
                if((SIGN == 1) && any(iplot < -istdn))
                    %     go = 0;
                end;
                if((SIGN == 0) && any(iplot > istdp))
                    %   go = 0;
                end;
                if(max(abs(iplot)) > par.t_maxamp)
                    fprintf(1, 'Candidate event rejected for being too large: %7.3f vs %7.3f max\n', ...
                        max(abs(iplot)), par.t_maxamp);
                    go = 0; % reject all the events that are larger than this...
                    iplot=[];
                    continue; % this also ends the analysis of this trace...
                end;
                if(go == 1) % if event is acceptable...
                    %                   fprintf(1, 'event %d accept\n', i);
                    if(SIGN == 1) % looking at positive-going events
                        [imax, jmax] = max(iplot);
                        m = m + 1;
                        ihist(m)  =  imax;
                        qhist(m) = sum(TRACE(k,qwin)-bl);
                        iamp = ihist;
                        i95 = decay_start*imax;
                        jmx = floor((length(iplot) - jmax)*0.9) + jmax;
                        [i95f p95f] = findnear(iplot(jmax:jmx), i95); % find point 95% way down falling slope
                        if(length(p95f) == 2)
                            p95f = p95f(2)+jmax - 1;
                        else
                            p95f = p95f + jmax - 1;
                        end;
                        if(fit_re_peak)
                            tpk = twin(jmax);
                        else
                            tpk = 0;
                        end;
                        time_fit = twin(p95f:jmx)-tpk;
                        if(jmx-p95f < 2)
                            tau = NaN;
                        else

                            [a0 a1 tau] = expfit(twin(p95f:jmx)-tpk, iplot(p95f:jmx)); % chebyshev fit.....
                            %         fprintf(1, 'expfit: %f +  %f *exp(-t/%f)\n', a0, a1, tau);
                            [fpar, chisq, niter, volt_fit1] = mrqfit('exponential', [a0 a1 tau], twin(p95f:jmx)-tpk, iplot(p95f:jmx), [], ...
                                [], [], [], 50, []);
                            %         fpar
                            if(isempty(volt_fit1))
                                fiti = a0 + a1*exp(-time_fit/tau);
                                fitmeth = 'c'; % chebyshev
                            else
                                a0 = fpar(1); a1=fpar(2); tau = fpar(3);
                                fiti = a1*exp(-(twin(p95f:jmx)-tpk)/tau)+a0;
                                fitmeth = 'm'; % marquardt
                            end;

                        end;
                        if(tau < p || tau > 2*iplus)
                            tau = NaN;
                        end; % limit taus to reasonable/acceptable values
                        itau(m) =  tau;
                        itaum(m) = 'm';
                        % fiti = a1*exp(-tn(p95f:end)/tau)+a0;
                    else % looking at negative going events
                        [imin, jmin] = min(iplot);
                        m = m + 1;
                        ihist(m) = imin;
                        iamp = ihist;
                        qhist(m) = sum(TRACE(k,qwin)-bl);
                        if(fit_re_peak)
                            tpk = twin(itemp_max);
                        else
                            tpk = 0;
                        end;
                        jmx = floor((length(iplot) - jmin)*0.9) + jmin;
                        time_fit = twin(p95f:jmx)-tpk;
                        if(jmx-p95f < 2)
                            tau = NaN;
                        else
                            %     [a0 a1 tau] = expfit(twin(p95f:jmx)-tpk, -iplot(p95f:jmx)); % chebyshev fit.....
                            %    fprintf(1, 'expfit(minus): %f +  %f *exp(-t/%f)\n', a0, a1, tau);
                            %    if(tau < 0)
                            %        tau = 10;
                            %    end;
                            a0 = mean(iplot(1:10));
                            if(ucmd == 2)
                                a1 = iscamp(i);
                            else
                                a1 = min(iplot);
                            end;
                            tau = par.t_decay;
                            %   a0
                            [fpar, chisq, niter, volt_fit1] = mrqfit('exponential', [a0 a1 tau], twin(p95f:jmx)-tpk, iplot(p95f:jmx), ...
                                [], [0,1,1],...
                                [-1000, -10000, 0.2], [1000, 0, 200], 500, []);
                            %      fpar

                            if(isempty(volt_fit1))
                                fiti = (a0 + a1*exp(-time_fit/tau));
                                fitmeth = 'c';
                            else
                                a0 = fpar(1); a1=fpar(2); tau = fpar(3);
                                fiti = a0 + a1*exp(-time_fit/tau);
                                fitmeth = 'm';
                            end;
                        end;
                        if(tau < p || tau > 2*iplus)
                            tau = NaN;
                        end; % limit taus to reasonable/acceptable values
                        itau(m) = tau;
                        itaum(m) = 'm';
                    end;

                    figure(hmini);
                    line(twin, iplot);
                    hold on;
                    if(fitmeth == 'c')
                        plot(time_fit, fiti, 'g');
                    else
                        plot(time_fit, fiti, 'r');
                    end;

                    if(m == 1)
                        icur = zeros(10000,length(iplot));
                    end;
                    icur(m,:) = iplot;
                    icstart(m) = icand(i);
                    icrec(m) = k;
                    icwin(m) = length(iwin);
                    ibl(m) = bl;

                    if(cmd == 2)
                        isamp(m) = iscamp(i);
                        icoff(m) = ioff(i);
                        ibl(m) = bl;
                    end;
                    hold on

                    evcount = evcount + 1; % keep track of events (even if they are not included in traces...)
                    if(evcount > 1)
                        isihist(m) = time(j)-time(lj);  % interevent distribution
                    end;
                    ieventtime{k}(m) = time(j);
                    if(exist('iplot', 'var') && ~isempty(iplot))
                        ieventamp{k}(m) = max(abs(sgolayfilt(iplot, 4, 11)));
                    else
                        ieventamp{k}(m) = NaN;
                    end;
                end; % of acceptable events (go = 1)
            else
                fprint(1, 'Candidate Event %d rejected for not fitting in window\n', j);
            end; % of window size check (does it fit?
        else
            fprintf(1, 'Rejecting candidate %d due to overlap n:%d p:%d\n', j, noverlap, poverlap);
        end; % of first level candidate event qualification check
        lj = j; % keep running track of most recent event
        jl = j+floor(iplus/p);
    end; % of the candidate event search loop (i)
end; % of the records loop (k)
if m == 0
    QueMessage('No Events detected', 1);
    analyzeflag = 0;
    return;
end
QueMessage(sprintf('EPSC mini analysis: Total = %d events identified',m));
icur = icur(1:m,1:end);

fprintf(1, 'EPSCMini - found %d events\n', m);
for i = 1:length(ieventamp)
    okevs = find(~isnan(ieventamp{i}));
    ieventtime{i} = ieventtime{i}(okevs);
    ieventamp{i} = ieventamp{i}(okevs);
    s = size(ieventamp{i});
    s(2)
end;

if(cmd == 2)
    tempw = [-predel:p:p*(length(tmpl)-floor(predel/p)-1)]+p;

    if(length(tempw) > length(tmpl))
        tempw = tempw(1:length(tmpl)); % just clip it to match.
    end;
end;
%if(par.t_disp)
if(0)    % plot the individual events... in multiplot format (all of them!)
    ng = ceil(sqrt(size(icur, 1)));
    hmp = multiplot('miniEPSCs', ng, ng, 1);
    mini = 0; maxi = 0;
    for i = 1:m
        axes(hmp(i));
        plot(twin, icur(i,:), '-k');
        if(cmd == 2)
            hold on
            plot(tempw, tmpl*isamp(i), '-r'); % scale...
            hold off
        end;
        mi = min(icur(i,:));
        if(mi < mini)
            mini = mi;
        end;
        ma = max(icur(i,:));
        if(ma > maxi)
            maxi = ma;
        end;
        title(sprintf('%d', i), 'fontsize', 7);
    end;
    mit = min(twin);
    mat = max(twin);
    for i = 1:m
        axes(hmp(i));
        set(gca, 'YLim', [mini maxi]);
        set(gca, 'XLim', [mit, mat]);
    end;
end;

% PSC_Scaled_meanvar(icur, isamp, m); % call sepately.
%
%--------------------------------------------------
%
h = newfigure('Mini_analysis', 'Mini Analysis Figure'); % check for pre-existing window
figure(h); clf;
set(gca, 'fontsize', fsize);
xlabel('ms');
ylabel('pA');
datac_setcrosshair(gca, 'Tr2_plot', 'ms', 'pA', grcxpos(1,:));

%subplot('Position', [0.19+2*0.8/3 0.12+0.8/3 0.8/3 0.8/3]);
%subplot(3,3,2);
subplot('Position', grpos(2,:));
cla;
icurn = icur(1:m,:);
icur = icurn;
sz = size(icur);

med = median(ihist); % find median value.
nbig = 0;
nsmall = 0;

for i = 1:sz(1)
    if(SIGN == 1)
        icur2(i,:)=icur(i,:)/max(icur(i,:));
    else
        icur2(i,:)=icur(i,:)/min(icur(i,:));
    end;

    if((SIGN == 1 && max(icur(i,:)) > med) || ...
            (SIGN == 0 && min(icur(i,:)) < med))
        nbig = nbig + 1;
        icurbig(nbig,:) = icur2(i,:); % upper half TRACEs
    else
        nsmall = nsmall + 1;
        icursmall(nsmall,:) = icur2(i,:); % lower half TRACEs
    end;
end;
if(sz(1) > 1)
    im = mean(icur, 1);
else
    im = icur;
end;
if(SIGN == 1)
    im = im/max(im);
else
    im = im/min(im);
end;
h=plot([twin(1) twin(end)], [0 0], '--');
set(h, 'color', [0.3 0.3 0.3]);
hold on;
if(SIGN == 0)
    h=plot([0 0], [-1.2 0.2],  '--');
    set(h, 'color', [0.3 0.3 0.3]);
else
    h = plot([0 0], [-0.2 1.2],  '--');
    set(h, 'color', [0.3 0.3 0.3]);
end;

if(nbig > 0)
    imbig = mean(icurbig);
    %   plot(twin, imbig, 'cyan');
end;
if(nsmall > 0)
    imsmall = mean(icursmall);
    %   plot(twin, imsmall, 'blue');
end;

% now analyze the psc shape. 20-80% rise time, half-width, and decay time constant.
% note that because we just analyze the mean trace, signs don't matter (its already normalized)
dt = mean(diff(twin)); % mean dt value
dtn = dt/100; % interpolation to finer time base
tn = [min(twin):dtn:max(twin)];
imn = interp1(twin, im, tn, 'spline'); % interpolate to finer time base
[imax, jmax] = max(imn); % find the max point.
i80 = 0.8*imax;
i20 = 0.2*imax;
[i80v p80] = findnear(imn(1:jmax), i80);
[i20v p20] = findnear(imn(1:jmax), i20);

rt = tn(p80)-tn(p20); % 20-80 rise time
rslope = (i80v-i20v)/rt; % and matching slope
i50 = 0.5*imax;
[i50r p50r] = findnear(imn(1:jmax), i50);
[i50f p50f] = findnear(imn(jmax:end), i50);

p50f=p50f+jmax-1;

thalf = tn(p50f)-tn(p50r); % width at half height.
i95 = decay_start*imax;
[i95f p95f] = findnear(imn(jmax:end), i95); % find point 95% way down falling slope

p95f = p95f+jmax - 1;

if(fit_re_peak)
    tpk = tn(jmax);
else
    tpk = 0;
end;
%[a0 a1 tau] = expfit(tn(p95f:end)-tpk, imn(p95f:end)); % chebyshev fit.....
it = mean(itau(find(~isnan(itau))));
mean(it)
[fpar, chisq, niter, volt_fit1] = mrqfit('exponential', [0 -1 it], tn(p95f:end)-tpk, imn(p95f:end), [], ...
    [], [], [], 50, []);
if(isempty(volt_fit1))
    fiti = a0 + a1*exp(-(time_fit-tpk)/tau);
    fitmeth = 'Chebyshev';
else
    a0 = fpar(1); a1=fpar(2); tau = fpar(3);
    fiti = a1*exp(-(tn(p95f:end)-tpk)/tau)+a0;
    fitmeth = 'Marquardt';
end;
plot([tn(p80) tn(p20)]-tpk, -[imn(p80) imn(p20)], 'o', 'markersize', 1.5, 'markerfacecolor', 'red', 'markeredgecolor', 'red');
plot(tn(p95f:end)-tpk, -fiti, '-r');
if(cmd == 2)
    plot(tempw-tpk, -tmpl/max(tmpl), '-b'); % draw template also
end;

plot(twin-tpk, -im, 'green'); % data on top of the template
plot(tn(p95f:end)-tpk, -fiti, '-r'); % put the fit ON TOP of the data

set(gca, 'XLim', [min(twin) max(twin)]);
set(gca, 'YLim', [-1.2 0.2]);
set(gca, 'fontsize', fsize);
ylabel('Normalized TRACE');
xlabel('ms');
datac_setcrosshair(gca, 'Tr2_norm', 'ms', 'I', grcxpos(2,:));

% histograms
%subplot('Position', [0.05 0.07 0.8/4 0.8/3]);
% Amplitude distribution
%subplot(3,3,4);
subplot('Position', grpos(4,:));
cla;
if(NBINS <= 4)
    NBINS = 32;
end;
[nh, xh] = hist(-ihist, NBINS); % 64 bins...
bar(xh, nh);
set(gca, 'fontsize', fsize);
xlabel('Amplitude (pA)');
ylabel('# events');
datac_setcrosshair(gca, 'Tr2_hist', 'pA', 'N', grcxpos(4,:));

m_i = mean(ihist);
s_i = std(ihist);

[qnh, qxh] = hist(-qhist, NBINS); % 64 bins...


% charge distibution
%subplot('Position', [0.10+0.8/4 0.07 0.8/4 0.8/3]);
%subplot(3,3,7);
subplot('Position', grpos(7,:));
cla;
bar(qxh, qnh);
set(gca, 'fontsize', fsize);
xlabel('Q (fC)');
ylabel('# events');
datac_setcrosshair(gca, 'Tr2_qhist', 'Q', 'N', grcxpos(7,:));

m_q = mean(qhist);
s_q = std(qhist);

[inh, ixh] = hist(isihist, par.nbins); % 64 bins...


% interevent interval
%subplot(3,3,5);
%subplot('Position', [0.15+2*0.8/4 0.07 0.8/4 0.8/3]);
subplot('Position', grpos(5,:));
cla;
bar(ixh, inh);
set(gca, 'fontsize', fsize);
xlabel('Interevent interval (ms)');
ylabel('# events');
datac_setcrosshair(gca, 'Tr2_ihist', 'ms', 'N', grcxpos(5,:));

% Amplitude versus decay tau
%subplot(3,3,8);
%subplot('Position', [0.20+3*0.8/4 0.07 0.8/4 0.8/3]);
subplot('Position', grpos(8,:));
cla;
io = find(itaum == 'm');
plot(itau(io), abs(iamp(io)), 'ro', 'markersize', 2.0);
hold on
ix = find(itaum == 'c');
plot(itau(ix), abs(iamp(ix)), 'gx', 'markersize', 2.0);
set(gca, 'fontsize', fsize);
ylabel('mEPSC (pA)');
xlabel('Decay Tau (ms)');
datac_setcrosshair(gca, 'Tr2_tauhist', 'ms', 'pA', grcxpos(8,:));

m_e = mean(isihist);
s_e = std(isihist);

% amplitude versus time, all events
%subplot(3,3,6);
subplot('Position', grpos(6,:));
cla;
et = [ieventtime{:}];
et = et(find(~isnan(et)));
ea = [ieventamp{:}];
ea = ea(find(~isnan(ea)));
plot(et, ea, 'ro', 'markersize', 2.0);
set(gca, 'fontsize', fsize);
ylabel('mEPSC (pA)');
xlabel('Time (ms)');
datac_setcrosshair(gca, 'Tr2_eventime', 'ms', 'pA', grcxpos(6,:));

% plot cumulative amplitude distribution...
%subplot(3,3,9);
subplot('Position', grpos(9,:));
cla;
ea = [ieventamp{:}]; % get all the event amplitudes.
xa = [0:2:50*ceil(max(ea)/50)];
eas = sort(ea);
na = cumsum(histc(eas, xa));


if(~isempty(na))
    stairs(xa, na, 'b-');
    set(gca, 'fontsize', fsize);
    ylabel('Cumulative events');
    xlabel('Amp (pA)');
    datac_setcrosshair(gca, 'Tr2_cumulative', 'ms', 'pA', grcxpos(9,:));
end;



EPSC.channel = par.channel; % save the whole parameter structure too.
EPSC.rise = par.t_rise;
EPSC.decay = par.t_decay;
EPSC.thresh = par.t_thresh;
EPSC.lpf = par.t_lpf;
EPSC.tstart = par.tstart;
EPSC.tend = par.tend;
EPSC.maxamp = par.t_maxamp;
EPSC.eventtime = ieventtime;
EPSC.eventamp = ieventamp;
EPSC.rt=rt;
EPSC.rslope=rslope;
EPSC.thalf = thalf;
EPSC.a0 = a0;
EPSC.a1 = a1;
EPSC.tau = tau;
EPSC.Imean = m_i;
EPSC.Istd = s_i;
EPSC.Qmean = m_q;
EPSC.Qstd = s_q;
EPSC.ISImean = m_e;
EPSC.ISIstd = s_e;
EPSC.tavg = twin;
EPSC.iavg = im;
EPSC.itau = itau;
EPSC.icur = icur;
EPSC.icstart = icstart;
EPSC.icwin = icwin;
EPSC.icrec = icrec;
EPSC.ibl = ibl;
EPSC.isamp = isamp;
EPSC.nevents = m;

switch(cmd)
    case 1
        EPSC.Method = 'Threshold';
    case 2
        EPSC.Method = 'Clements&Bekkers';
    otherwise
end;
if(fit_re_peak == 1 || cmd == 1)
    EPSC.Align='peak';
else
    EPSC.Align='foot';
end;

dcmd = sprintf('CONTROL(sf).EPSC%d=EPSC;', par.channel);
QueMessage(sprintf('Data stored in control as EPSC%d', par.channel), 1);
eval(dcmd);
CONTROL(sf).Commentary = sprintf('EPSC_mini: %.1f-%.1f ms', par.tstart, par.tend);

% text box...
%subplot('Position',[0.14+2*0.8/3, 0.13+2*0.8/3, 0.2, 0.35])
subplot('Position', grpos(3,:));
%subplot(3,3,3);
cla;
axis([0,1,0,1])
axis('off')
ht(1)=text(0,0.9,sprintf('%-12s R[%d:%d]',DFILE.filename, DFILE.frec, DFILE.lrec));
ht(2)=text(0,0.82,sprintf('gain:%4.1f  LPF:%4.1f kHz', DFILE.igain, DFILE.low_pass(1)));
ht(3)=text(0,0.74,sprintf('Protocol: %s', CONTROL(sf).protocol));
set(ht(3), 'interpreter', 'none'); % un-TeX the line - this is a filename and won't have tex chars, but might have an underscore.
ht(4)=text(0,0.66,sprintf('I:   \\mu=%5.1f \\sigma=%5.1f pA CV=%5.2f (N=%d)', m_i, s_i, s_i/m_i, length(ihist)));
ht(5)=text(0,0.58,sprintf('Q:   \\mu=%5.1f \\sigma=%5.1f pA CV=%5.2f (N=%d)', m_q, s_q, s_q/m_q, length(qhist)));
ht(6)=text(0,0.50,sprintf('IEI: \\mu=%5.1f \\sigma=%5.1f pA CV=%5.2f (N=%d)', m_e, s_e, s_e/m_e, length(isihist)));
ht(7)=text(0,0.42,sprintf('Decay \\tau=%6.3f ms  20-80RT = %6.3f ms  ', tau, rt));
ht(8)=text(0,0.34,sprintf('Fit Method: %s  Align: %s Detect: %s', fitmeth, EPSC.Align, EPSC.Method));
ht(9)=text(0,0.28,sprintf('Half-width=%6.3f ms  ', thalf));
atau = itau(find(~isnan(itau))); %#ok<FNDSB>
ht(10)=text(0,0.20,sprintf('\\tau:  \\mu=%5.3f ms \\sigma=%5.3f ms N=%4d', mean(atau), std(atau), length(atau)));
ht(11) = text(0, 0.12, sprintf('tstart: %7.1ms  tend: %7.1ms', par.tstart, par.tend));
ht(12)=text(0.75,0.03,sprintf('%s',date));
for i = 1:length(ht)
    set(ht, 'Fontsize', 9);
    %set(ht, 'Fontname', 'helvetica');
end;
orient landscape
analyzeflag = 0;
return;


